#include "EDO.h"
#include <iostream>

EDO1::EDO1(int a)
{
N=a;
}


double EDO1::rk2 (double x0,double y0,double xf)
{
double f0=evalyp(x0,y0);
double h = (xf-x0)/N;
double k0=h*f0;
double xi=x0;
double yi=y0;

for(int i=0; i<N; i++)
{
double fi=evalyp(xi,yi);
double ki=h*fi;
yi=yi+h*evalyp(xi+h/2,yi+ki/2);
xi=xi+h;
}
return yi;
}


double EDO1::rk4(double x0,double y0,double xf)
{
double h=(xf-x0)/N;
double xi=x0;
double yi=y0;

for(int i=0; i<N; i++)
{
double fi=evalyp(xi,yi);
double k1=h*fi;
double k2=h*evalyp(xi+h/2,yi+k1/2);
double k3=h*evalyp(xi+h/2,yi+k2/2);
double k4=h*evalyp(xi+h, yi+k3);

yi=yi+(1.0/6)*(k1+2*k2+2*k3+k4);
xi=xi+h;
}
return yi;
}



double EDO1systeme::getyf()
{
return yf;
}


double EDO1systeme::getzf()
{
return zf;
}



void EDO1systeme::rk4(double x0, double y0, double z0, double xf)
{

double h=(xf-x0)/N;
double xi=x0;
double yi=y0;
double zi=z0;
for(int i=0; i<N; i++)
{
double fi= evalyp(xi,yi,zi);
double gi= evalzp(xi,yi,zi);
double k1= h*fi;
double l1= h*gi;
double k2= h*evalyp(xi+h/2,yi+k1/2,zi+l1/2);
double l2= h*evalzp (xi+h/2,yi+k1/2,zi+l1/2);
double k3= h*evalyp(xi+h/2,yi+k2/2,zi+l2/2);
double l3= h*evalzp(xi+h/2,yi+k2/2,zi+l2/2);
double k4= h*evalyp(xi+h,yi+k3,zi+l3);
double l4= h*evalzp(xi+h,yi+k3,zi+l3);
yi=yi+(1.0/6)*(k1+2*k2+2*k3+k4);
zi=zi+(1.0/6)*(l1+2*l2+2*l3+l4);
xtab=xi;
ytab=yi;
yptab= evalyp(xi,yi,zi);
ztab= zi;
zptab= evalzp(xi,yi,zi);
xi=xi+h;
}
yf=yi;
zf=zi;
}



double EDO2:: getypf()
{
return yf;
}


double EDO2::evalyp (double x, double y,double z)
{
return z;
}


double EDO2::evalzp (double x, double y,double z)
{
return evalys(x,y,z);
}




EDO1systeme::EDO1systeme (int n)
{
n=N;
xtab= new double[n];
ytab= new double[n];
yptab= new double [n];
ztab = new double[n];
zptab= new double [n];
}

EDO1systeme::~EDO1systeme()
{
delete []xtab;
delete []ytab;
delete []yptab;
delete []ztab;
delete []zptab;
}

void EDO1systeme::ecrirexyz(const char *tp2)
{
ofstream FichierSortie;
FichierSortie.open(tp2);
FichierSortie << setw(9)<<xtab<< setw(20)<<ytab<<setw(20)<<yptab<< setw(20)<<ztab<<setw(20)<<zptab<<endl;
}

# include <iostream>

using namespace std;

class EDO1
{
public:
int N;
EDO1(int N);
virtual double evalyp(double x,double y)= 0;
double rk2(double x0,double y0, double xf);
double rk4(double x0, double y0, double xf);
};

class EDO1systeme
{
public:
double yf;
double zf;
int N;
double*xtab;
double*ytab;
double*yptab;
double*ztab;
double*zptab;


EDO1systeme(int n);
~EDO1systeme();
double getyf();
double getzf();
virtual double evalyp (double x, double y,double z)=0;
virtual double evalzp (double x, double y,double z)=0;
void rk4(double x0, double y0, double z0, double xf);
};

class EDO2: public EDO1systeme
{
public:
EDO2(int n):EDO1systeme(n){};
double evalyp (double x, double y,double z);
double evalzp (double x, double y,double z);
virtual double evalys(double x, double y, double z)=0;
double getyf();
double getypf();
};

# include <iostream>
# include "EDOtest.h"


double EDO1test::evalyp(double x,double y)
{
double f=-y;
return f;
}

/double EDO1systemetest::evalyp(double x,double y,double z)
{
double f=z;
return f;
}

double EDO1systemetest::evalzp(double x,double y,double z)
{
double g=-y;
return g;
}


double EDO2test::evalys(double x, double y, double z)
{
double h= -y;
return h;
}

# include <iostream>
# include "EDO.h"

using namespace std;


class EDO1test:public EDO1
{
public:
EDO1test(int n): EDO1(n) {};
double evalyp(double x,double y);
};

class EDO1systemetest: public EDO1systeme

{
public:
EDO1systemetest(int n): EDO1systeme(n) {};
double evalyp (double x, double y,double z);
double evalzp (double x, double y,double z);
};




class EDO2test:public EDO2
{
public:
EDO2test(int n): EDO2(n) {};
double evalys(double x, double y, double z);
};

# include <iostream>
# include <iomanip>
# include "EDOtest.h"
# include <cmath>

using namespace std;

int main()
{

EDO1test A(10);
EDO1test B(100);
EDO1test C(1000);
cout<<setprecision(12) <<A.rk2(0,1,2)<<endl;
cout<<setprecision(12) <<B.rk2(0,1,2)<<endl;
cout<<setprecision(12) <<C.rk2(0,1,2)<<endl;
cout<<setprecision(12) <<A.rk4(0,1,2)<<endl;
cout<<setprecision(12) <<B.rk4(0,1,2)<<endl;
cout<<setprecision(12) <<C.rk4(0,1,2)<<endl;

EDO1systemetest D(100);
EDO1systemetest E(1000);
cout<<D.rk4(0,0,0,5*M_PI/4.0)<<endl;
cout<<E.rk4(0,0,0,5*M_PI/4.0)<<endl;

EDO2test F(100);
cout<<F.rk4(0,0,0,5*M_PI/4.0)<<endl;

return 0;
}